Competition between Diffusion and Fragmentation: 
An Important Evolutionary Process of Nature 
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We investigate systems of nature where the common physical processes diffusion and fragmen- 
tation compete. We derive a rate equation for the size distribution of fragments. The equation 
leads to a third order differential equation which we solve exactly in terms of Bessel functions. The 
stationary state is a universal Bessel distribution described by one parameter, which fits perfectly 
experimental data from two very different system of nature, namely the distribution of ice crystal 
sizes from the Greenland ice sheet and the length distribution of Q-helices in proteins. 



A diffusion process is one of the most important 
and common physical phenomena of nature. Coherent 
structures, like crystals and structural elements of bio- 
molecules may shrink or grow gradually and randomly in 
a way that resembles a diffusive motion of their bound- 
aries. On the one hand, inhomogeneities in a particular 
structure may disappear or be resolved into neighboring 
structures due to diffusion. On the other hand, diffu- 
sion typically tends to make the largest of the coherent 
structures larger. In contrast to this gradual modifica- 
tion of individual fragments, one might in addition en- 
counter a completely different and abrupt physical phe- 
nomena, namely that of fragmentation. Fragmentation 
occurs in for example growing ice crystals, when the dif- 
fusively growing crystal is subjected to stresses that can 
cause a piece to break off, thus leading to two individual 
crystals, each of them continuing the competitive process 
between diffusion and fragmentation. 

We believe that such interplay between gradual diffu- 
sion and rare but drastic fragmentation is a very common 
phenomenon in nature. In this letter we derive a general 
rate equation for the dynamical evolution of the density 
of fragments of a given size. The rate equation consists of 
a local diffusion term and non-local fragmentation terms. 
Differentiating once we obtain a third order differential 
equation in the size distribution of fragments. We solve 
the equation exactly by means of an eigenvalue expansion 
in terms of Bessel functions. The different eigenvalues cor- 
respond to different decay times such that in the infinite 
time limit we obtain the stationary distribution, corre- 
sponding to zero eigenvalue. This distribution is char- 
acterized by one parameter (the ratio between the diffu- 
sion and the fragmentation constants) and constitutes a 
universal stationary distribution for these competing pro- 
cesses. We apply experimental data from two very differ- 
ent physical processes, the size distribution of ice crystals 
in the Greenland ice sheet Q] and the length distribution 
of a-helices in proteins of low homology. In both cases we 
find that the universal distribution fits the experimental 
data perfectly. In addition, for the ice crystals, we are 



able to map out the entire dynamics of ice crystals 2000 
years back in time 0. 

In order to clarify the description, we consider a one- 
dimensional simplification of the process in which the sizes 
of for instance ice-crystal grains are projected onto the co- 
ordinate x. The mechanism we propose then is described 
in terms of the density N(x, t) of "objects" (e.g. ice crys- 
tals or a-helices in proteins) of a given length x at a given 
time t. At a given time it is possible to increase or de- 
crease N(x,t) by diffusion characterized by the diffusion 
constant D. Also, N(x,t) can be altered by fragmen- 
tation of the "objects" characterized by a fragmentation 
constant /, defined as the average number of break-ups 
in a time interval dt over a length L. Thus, the number 
of fragments is given by fLdt. The constants D and / 
can depend on temperature and on quantities which are 
characteristics of the "objects" . The resulting equation is 

dN(x,t) _ D d 2 N{x,t) 



dt 



dx 2 



-fxN(x,t) + 2f / dx'N(x,t). 

J X 

(1) 

Note that by this equation, the total "mass" of the frag- 
ments J Q xN{x, t)dx is conserved at any time in the evo- 
lution. The first term on the right hand side represents 
diffusion, the second term gives the fragmentation of ob- 
jects with length x, and the last term is the contribution 
from fragmentation of objects larger than x. The factor 
of two in the last term follows formally from the condition 
of mass conservation, and results from the fact that there 
are two distinct ways to cut an object of length x + x' into 
segments of length x and x' . Similar approach has been 
applied in 2] for the stochastic dynamics of microtubules 

a 

We introduce the dimensionless variables x — x /xq and 
t = t/to, where 



t = (Dfr^ and xo = (D/fy/ 3 



(2) 
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This amounts to considering eq. QJ with D = f = 1, 
which we shall do in the following. We now differentiate 
eq. Q with respect to x and obtain 

d t d x N(x, t) = 8%N(x, t) - W(x, t) - xd x N(x, t). (3) 

This equation can be solved by separation, using the 
boundary condition that N(x, i) — ► for x — > oo. The 
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result can be expressed in terms of the Airy function 

f°° ( k 3 ' 

J dk cos ykx + — 



A(x) 



(4) 



which in turn can be expressed in terms of Bessel func- 
tions [4|, 



Ji 



2a; 3 / 2 
3 

2b| 3 / 2 



for x > 
7 ,'2bl 3 / 2 



(5) 



for x < 0. 
The solution of © is then 

N(x, t)=J2 C n e x " f B{x + A„), B{x) = -d 2 A(x), (6) 

where the sum may be an integral in case the eigenvalues 
A n are in a continuous range. The function B(x) is thus 
related to the Airy function by the second eq. ©. One 
further has 



B{x) = -xA(x). 



(7) 



This follows from the second eq. 10 by performing two 
differentiations of the expressions © and by use of stan- 
dard Bessel function recursion relations. For a discussion 
of the Airy function we refer to Q. 




FIG. 1: An example showing the time development of 
N(x,t) with a secondary peak, which ultimately disappears. 
The specific values of the expansion parameters C„ are 
(C ,C\,C 2 ...) = (-7, -.2, -.2, .05, -.035) . 

We now impose the boundary condition that there 
should be no objects of size x = 0, 



N(0,t) = 0. 



(8) 



To implement the requirement JSJ) on the solution © we 
thus need that the eigenvalues A„ are the zeros of B(x). 
From eq. JJJ we see that x = is a zero. For x > there 
are no zeros, whereas for x < there exist an infinite 
number of zeros given as the solution of the equation 



Ji 



2 1 A„ 1 3/2 



.]_ 



2|A„| 3 / 2 



0. 



(9) 



These zeros can easily be found numerically. The first few 
are given approximately by Ao = 0, Ai = —2.338, A2 = 
-4.088, A 3 = -5.521, A 4 = -6.787, A 5 = -7.945, A 6 = 
-9.023. 

In order to find the the constants C n in the solution @) 
we notice that the basic equation © is of third order, and 
hence does not lead to orthogonal functions. Instead we 
notice that from J7J) the Airy function satisfies the second 
order differential equation, 



%A(x) = xA(x), 



(10) 



for which standard methods are applicable: We define the 
auxiliary function 



M{x,t) = J2 c ne Xnt A(x + \ n ), 



(11) 



n=0 



and we can then easily show from the second order dif- 
ferential equation i|lU[l that the functions A(x + A„) are 
orthogonal, and hence C n can be obtained in terms of an 
integral over the initial function M(x, 0) times A(x + A„). 
Since, however, we also have 



d 2 x M(x,t) = -N(x,t) 



(12) 



from the initial function M(x,0) can be expressed in 
terms of the corresponding initial function N(x, 0), giving 
the result 



C - i 



dx A(x + A n ) 



dx' {x - x')N{x', 0) - C A(x) 



(13) 



where 



and Cn = 



dxA(x 
-1 



A 



A(Q) 



,f = A\\ n )\ 
dxxN(x,0). 



(14) 



This then allows us to write down the solution N(x, t) for 
any given initial function N(x, 0). A characteristic feature 
of this solution is that the mean value (x) saturates for 
large times. Using the solution one finds 



(x) 



J™xN{ Xl t) 3 2 / 3 r(4/3) 



f~N(x,t) t- 



r(2/3) 



(15) 



which means that (x) — > 2.096(Z3//) 1 / 3 . We also find 
that the dispersion saturates, (x 2 — (x) 2 ) — » 0.6074(x) 2 , in 
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marked contrast to pure diffusion, where the dispersion in- 
creases without bound. The time development of N(x, t) 
is governed by the Boltzmann-like factors exp(A„<). In 
Fig. 1 we have shown by an example that these factors 
can give rise to additional secondary bumps, which ulti- 
mately disappear for larger times. It would be interesting 
to find actual examples where these structures are ob- 
served experimentally. 
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FIG. 2: The black dots represent the length distribution of 
a-helices from 299 high resolution protein structures with low 
homology [5j. The error bars are estimated as the square root 
of the numbers. The full curve is the stationary distribution 
(i.e. infinite time limit) obtained from eq. Q plotted versus 
the helix length measured in terms the number of amino acids. 
The ratio between the diffusion constant D and the fragmenta- 
tion rate / resulting in the best fit to the data is (y) 1 ^ 3 ~ 6.1 
number of amino acids. Note that N(x) vanishes for a helix 
length equal to 3, as a a-helix needs one turn to be identified, 
see Q. 



We now describe two widely different application of the 
suggested interplay between gradual diffusion and sudden 
fragmentation. The examples are selected from the con- 
straint imposed not so much by the generic nature of the 
two processes, but also by the requirement that random 
merging of structures should be insignificant. Further, 
the two examples are taken from cases where it is rea- 
sonable to assume that fragmentation occurs uniformly 
on the size axis x, as implicitly assumed in the basic eq. 
||TJ of the process. The two experimental data sets be- 
hind our studies are respectively the length distribution 
of a-helices in proteins, and the size distribution of ice 
crystals in ice sheets. In the former case we will be able 
to extract the ratio y by performing a least squares fit of 
the stationary state solution and in the latter case we do 
an exponential fit of (x) (t) using the two leading terms in 



the solution, 



D and / from J2J and i|15[l and by noting that the char- 
acteristic time t is related to the fragmentation by the 
second largest eigenvalue Ai, we obtain r = -l/(Ai/) Q. 
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FIG. 3: Distributions of ice crystal sizes at depths 115m, 165m, 
220m, 330m, 440m and 550m. The crystal size is defined as 
the vertical extension of the individual crystals. The ragged 
curves are the measured histograms and the smooth lines are 
the temporal evolution predicted by eq. Q starting from the 
initial distribution at 115m. The total counts of ice crystals 
decreases with depth (due to the overall increase of sizes) un- 
til the steady state is reached. The values of the diffusion 
constant and the fragmentation rate used in these plots are 
D = 2.8 • 10~ 3 mm 2 /year and / = 3.6 ■ 10~ 4 /cm ■ year. Note 
that the distributions vanish at a small, finite value to account 
for the experimental bias; see Q. 

The length distribution of a-helices [j| is taken from a 
database of 299 high resolution structures with low ho- 
mology extracted from the PDB (Protein Data Bank) p| . 
Secondary structure have been calculated on the basis 
of backbone hydrogen bonding and in total 2101 he- 
lices are used 0. Helices of lengths between 5 and 31 
amino acids have been counted resulting in the distribu- 
tion shown as circles in Fig. |3 Clearly, the distribution 
exhibits a maximum around a length of 7 amino acids 
followed by a long tail. The data do not fit well any sim- 
ple statistical distribution and a polynomial fit requires 
at least a fourth order expression to become acceptable 

5]. We consider this distribution as the result of an 'in- 
finitely' long evolutionary process where diffusive growth 
competes with intermittent fragmentation. The station- 
ary solution to eq. (JIJ is defined by one free parameter 
(the ratio D/f) and it nicely reproduces the overall rise 
and decline of the observed helix frequency with length 

7]. There are several points to make in this connection: 
1) Considering ensembles of a-helices, one effectively ran- 
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domizes the particular evolutionary advantages of growing 
or shrinking any particular a-helix in any particular pro- 
tein. Thus whereas a particular change may have a pur- 
pose, then the overall process is probably well represented 
by random domain adjustments and occasional fragmen- 
tation. 2) Though D and / have essentially the same 
molecular origin, simple energy arguments suggest that 
a point mutation in the middle of the helix to a residue 
not forming hydrogen bonds has a statistical weight of 
~ 1/400 relative to the weight of growing/shrinking the 
helix 0. This corresponds to the ratio (D//) 1 / 3 rs 7.4 in 
reasonable agreement with our result. 3) That in princi- 
ple other processes also acts to limit the maximum length 
of a-helices, as f.ex. the total size of proteins limits the 
maximal length of their sub elements. This maximum 
length will only influence the distribution of the very long 
a-helices. Overall we see the resulting distribution of a- 
helices as an ensemble average, revealing a process that 
leaves a statistical stationary distribution of structural el- 
ements in much the same way as binding energies between 
amino acids tend to influence the frequencies of neighbor- 
ing relations between them in an ensemble of proteins . 

For the ice crystal distribution, we apply data from the 
North Greenland Ice Core Project (NorthGPJP) which 
provides paleoclimatic information back to at least 115 
kyr before present (B.P.) [T^j. Each year, precipitation 
on the ice sheet covers it with a new layer of snow, which 
gradually transforms into ice crystals as the layer sinks 
into the ice sheet. The size distribution of ice crystals 
has been measured at selected depths in the upper 880 
m of the NorthGFJP ice core [ll|, which cover a time 
span of 5300 years pf . Fig. shows size distributions 



of ice crystals at selected depths down to 550 m. The 
crystal size x is defined as the vertical extent of a crystal 
(13| . Applying instead the horizontal extent of the crystal 
yields the same distribution curve thus supporting our 
assumption of using the linear size in our formalism. 

Each distribution exhibits a pronounced peak, indicat- 
ing a typical crystal size at each depth, followed by a 
decaying tail of relatively large crystals. The mean size 
becomes larger with depth and thus time until it saturates 
[l4iri5| . The distributions gradually change with time to- 
ward a universal curve, indicating a common underlying 
physical process in the formation of crystals. We have 
identified this process as an interplay between the frag- 
mentation of the crystals and the diffusion of their grain 
boundaries and is thus described by our general frame- 
work. 

In this Letter we have presented a general evolutionary 
scenario for the dynamics of objects subjected to a compe- 
tition between diffusion and fragmentation. The process 
results in an evolution equation of the density N(x,t) of 
objects of linear size x. In the stationary limit, this equa- 
tion predicts a "universal" distribution curve N(x) deter- 
mined only by the ratio y. This curve can be described 
exactly in terms of Bessel function leading to a power 
law uprise for small x followed by a tail of a stretched 
exponential form for large x (exp(~x^)). We speculate 
that this general competing principle might be the rele- 
vant mechanism in many other systems of nature and may 
result in distribution functions of a similar type. 

We are grateful to Mogens Levinsen and Anders Svens- 
son for discussions. 
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